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VERY RESTRICTED FOUR-BODY PROBLEM 

by 

Su-Shu Huang 
SUMMARY 

First, a state of motion of three finite bodies m l9 m 2 , m 3 
is idealized by an approximation to the law of mechanics 
such that m 2 and m 3 revolve around each other in circular 
orbits and that their center of mass revolves around m 1 
also in a circular orbit. The motion of a fourth body of an 
infinitesimal mass is then studied in a similar manner, as 
in the restricted three-body problem. 
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VERY RESTRICTED FOUR-BODY PROBLEM 
INTRODUCTION 

In another paper* the general behavior of an artificial satellite in the earth-moon- 
sun system was studied in terms of two three-body problems. In the present paper some 
justification will be provided for that approach by treating dynamically an idealized case 
of motion of an infinitesimal body of mass m in a system of three bodies mj , m 2 , and m 3 
so arranged that the center of mass, O', of m 2 and m 3 is revolving around the center of 
mass, 0, of the entire system in a circular orbit and m 2 and m 3 themselves are revolv- 
ing around 0' also in circular orbits. Such a state of motion of the three bodies is obvi- 
ously possible only in the form of approximation. However, if 

m l >> m 2 + m 3 >> m 

and if the separation A between and O' is very much greater than that separation a 
between m 2 and m 3 , both of these two conditions being true in the case of an artificial 
satellite in the earth-moon-sun system, the approximation will deviate from 

the actual solution of mechanics very little. 


AN INTEGRAL OF THE EQUATION OF MOTION 
FOR THE FOURTH BODY 


Further assume that the three bodies m u m 2 , m 3 always remain in the same plane 
and let the distances of mj and 0' from 0 be A x and A 2 , and those of m 2 and m 3 from 
0 ' be a x and a 2 . Now choose a rectangular coordinate system with its origin at 0 and 
its three axes ?, v , £ fixed in space, the $ -axis being perpendicular to the plane of the 
three finite bodies. Hence the coordinates of the four bodies may be written as 
m i(s r i> 7 7i>0) > .0), m 3 (<f 3,773,0), and m($, 77, £) . The equations of motion of the 

infinitesimal body m are given by 


dt 2 


~ Gm^ 






(i) 


*Huang, S.-S., "Some Dynamical Properties of the Natural and Artificial Satellites, 11 
NASA Technical Note D-50Z. 
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d 2 ?] 

dt 2 


v-Vx 

- Gmj —j- 
r i 


V “ 7 ? 2 

Gm 9 — “ Gm 3 — 

2 r/ r 3 3 


( 2 ) 


d 2 £ 

dt 2 




(3) 


where ri , r 2 , r 3 represent the distances of the infinitesimal body m from m x , m 2 , and 
m 3 , respectively. 


Let n x be the angular velocity with which the m 2 -m 3 system revolves around 0, and 
that with which m , revolves around O'. It can be easily seen that 


Cl = “ Ai cos Hjt , 

^2 = A 2 cos f 2 1 t 

7)2 = A 2 sin Qjt 

and 

C 3 - A 2 cos n x t 

t ] 3 = A 2 sin Qjt 


Tj ^ = - Aj^ sin Ojt ; (4) 

- a! cos fi 2 t » 

- aj sin 0 2 t ; (5) 


+ a 2 cos 0 2 t , 

+ a 2 sin Q 2 t . ( 6 ) 


To the approximation involved in the assumption of circular motions of the three 
finite bodies, we have, from the result of the two-body problem, 


and 


Obviously, 


<>1 


i_ 

G(mi +m 2 + ni 3 )”| 2 
A 3 


n 2 



A = A! + A 2 , and a = a! + a 2 . 


(7) 


( 8 ) 


(9) 


Next, choose anew rectangular coordinate system xyz centered at O' with the x- 
axis revolving with m 2 and m 3 and with the z-axis parallel to the 5-axis. The equations 
of transformation from the old one to the new one can easily be found to be 


C - A 2 cos fijt + x cos G 2 t - y sin H 2 t , 
V - A 2 sin Gjt + x sin Q 2 t + y cos 0 2 t , 

t = z • 


( 10 ) 

(11) 

( 12 ) 
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The coordinates of m 1} m 2 , and m 3 in the new system are given by 


and 


Xj = - A cos (Vl 2 - Qj) t , 

yj = A sin (ft 2 - flj) t , 

2 1=0; 

(13) 

X 2 = -a x , 

y2 - 0 . 

z 2 = 0 ; 

(14) 

x 3 = a 2 . 

Y3 = 0 . 

z 3 = ® ■ 

(15) 


Substituting Equations 4 through 6 and 10 through 12 in Equations 1 to 3, and utiliz- 
ing Equations 9 and 13 through 15 give the following result: 


d^x __ 

2n ^ : 
2 dt 

dt 2 

^ + 
dt 2 


d^z 

^ = 

■°”i 73 
F 1 

where 



a 2 n? > 
A a? x * 


A/) O-' 


Gm 


- xj x + a x x - a 2 

i —r ~ Gm 2 — — - Gm 3 — — - 

r l r 2 r 3 


l y "T^ yi r Gmi T? _Gm2 ^ 


3 ’ 


— Gm 


’ 


(16) 

(17) 

(18) 


(x-xj) 2 + (y-y x ) 2 + z 2 , 


(x + a^ 2 + y 2 + z 2 , 


r, = (x - a 2 ) + y* + z 


(19) 


If a function u is defined as 


U(x,y, z) S n. 


2 (x 2 + y*) 


2 '\ - 


a 2 n, 

t ^2 (x i x + yi y) 

u 2 


Gm 5 Gm 

+ 1 + — 

fi r 


2 0,11 3 


( 20 ) 


then Equations 16 through 18 assume the simplified form 


d 2 x dy BU 

dt* " 2n2 dT = ’ 


d 2 y dx _ BU 

+ 2n 2 Ht = ’ 


(21) 

( 22 ) 


d 2 z _ BU 


dt 


2 3z ’ 


(23) 
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which can be integrated to give, for each epoch of an infinitesimal time-interval, 

v 2 = 2U + constant , (24) 

where V is the magnitude of velocity in the xyz system of reference. Equation 24 plays 
a role in the present problem, just as Jacobi's integral in the restricted three-body 
problem. 


ZERO-VELOCITY SURFACES 

It follows from Equations 20 and 24 that the zero-velocity surface can be defined by 


“2 9 2 

(x 2 + y 2 ) 


A 2 9 Gm j Gm 2 

t n ! (x * x + y i y) + T7 + t: 


Gm q 


constant . 


(25) 


Since z x = o, 


x i x + yiy ~ rA cos ® » (26) 

where 0 is the angle between 0'm and 0' mi . In the case of the earth-moon-sun system, 
it is the angle subtended by the artificial satellite and the sun at O' . Thus, the zero- 
velocity surfaces are not fixed even in the rotating coordinate system; rather, they 
change with the position of m 1 . However, an instantaneous (or osculating) zero-velocity 
surface can be defined for each position of m x . It is in this sense that zero-velocity 
surfaces will be discussed. Indeed, the general behavior of the motion in the very 
restricted four-body problem can be understood by these osculating zero-velocity sur- 
faces just as that of the motion in the restricted three-body problem by the zero-velocity 
surfaces themselves. 

With the aid of Equations 7, 8, and 26, Equation 25 becomes 


1 

2 


(m 2 + m 3 )(x 2 + 


y 2 ) 


mjr 

A 2 


m j m 2 m 3 

cos © + — + + 

r l r 2 r 3 


constant . 


(27) 


Since the motion of m is of interest only when r « A, the term l/r x in Equation 27 can 
be expanded in terms of spherical harmonics, F n (cos 0); that is, 


r i 



P n (cos 0) 


(28) 


Taking only the first three terms in the right side of Equation 28 and substituting them in 
the place of l/r 1 in Equation 27, we obtain 
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(m 2 + m 3 )(x 2 + y 2 ) mjr 2 2m 2 2m 3 

+ 5 — (3 cos 0 - 1) + + = constant , (29) 

a 3 A 3 r 2 r 3 V ' 

where the term 2114 /A has been absorbed in the constant term. 

If a is now taken as the unit of length and m 2 + m 3 as the unit of mass, Equation 29 
reduces to 


(3 cos 2 ® 


1 ) 


+ 2(1 -M) 

r 2 



(30) 


where 


M 



(31) 


and c is a constant of integration. This differs from the zero-velocity surfaces of the 
restricted three-body problem only by the addition of a small perturbing term that con- 
tains the factor mj/A 3 . 


DOUBLE POINTS OF THE SURFACES 


Consider the change in position of the three double points l 19 l 2 , l 3 which are 
located on the x-axis when the perturbing term vanishes. Since this is now limited to 
the xy plane, 


0 = 6 - 6 0 , (32) 

where 0 and 0 Q are the respective angles that the positive x-axis makes with the vectors 
0 'm and O'mj. In the case of the earth-moon-sun system, 0 0 changes from 0 to in a 
period of the lunar month. Therefore, in a time-scale of a few hours 8 0 may be regarded 
as constant. 

Substituting Equation 32 into Equation 30, we obtain after reduction 

F(x,y) s (1 + /3)(x 2 +y 2 ) + 3/3 |~(x 2 -y 2 ) cos 2^ 0 + 2xy sin 2^ 0 1 


where 


2(1 -M) + 2m 
r 2 r 3 


C 


0 


(33) 


( 34 ) 


The conditions for double points are 
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\ |f- = [l + (1 + 3 cos 2 & 0 ) /3] x + 3</3 sin 2& 0 ) y 
(1 - “ * 2 > M(x - x 3 ) 


and 

5 Ip- = [l + (1 - 3 cos 2 & 0 ) p]y + 3(/3 sin 2e o ) x 

(i - m) y My _ „ 

r 3 ' r 3 

2 3 

from which all five double points can be determined on the xy plane. Double points are 
no longer the particular solution of the problem because they change with 6 0 and also 
because the problem is being treated only approximately (by taking only the first few 
terms in the series expansion of l/r 1} etc.). 

Since p is small, its second and higher orders can be neglected. It appears from 
Equation 36 that the .three double points which approach the x-axis when P = o have their 
y -coordinates of the order of P sin 2& 0 when P ^ 0. Thus, the term 3(/3 sin 2<9 0 )y in 
Equation 35 is of the order of (P sin 20 o ) 2 and can be neglected. Hence, Equation 35 is 
reduced to 


(35) 


(36) 


1 + (1 +3 cos 


) P : 


(1 - fi)(x - x 2 ) m(x - x 3 ) 


= 0 , 


(37) 


which differs from its counterpart in the restricted three-body problem only by the factor 

(1+3 cos 2& 0 ) P. 

Once the x -coordinates of the three double points are derived from the solution of 
Equation 37, their y-coordinates can be obtained by 

{%) M - F i sin 2e » ■ ( 38 > 

where 


which follows directly from Equation 36 except that now the values of x, r 2 , and r 3 are 
taken at one of the three double points Lj(i = 1,2,3) for the case p - o. 


The change in position of the three double points with P and 6 0 can be most con- 
veniently seen by first taking the derivatives of their coordinates with respect to p and 
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setting ,6 = 0. Consider the three points separately: (1) L 2 between + 00 and x 3 , (2) L x 
between x 3 and x 2 , and (3) l 3 between x 2 and -®. 


(1) Let the distance from m 3 in the x-direction to the double point L 2 be repre- 
sented by p. Then Equation 35 becomes, by neglecting the second and higher orders of 

P sin 2<? 0 , 

I" 1 + ( 1 + 3 cos 28 0 ) /fil ( 1 - p + p) - — — 7 ~ ~2 ~ -^2 - 0 ■ ( 39 ) 

L J (1 + P) P 

Differentiating Equation 39 with respect to p and setting p = o give 

0 = >2d + 3 COS 26 0^ • 

where 

(1 - m + p 0 )(1 + p 0 ) p 0 3 
2 3P 0 3 (1 + P 0 ) + 2p(l - P^) 

The symbol p 0 at the right side of Equation 41 is the solution of p for Equation 39 with 
p = o. Similarly the change in value of c which corresponds to the variation in position 
of L 2 can be computed. Differentiating Equation 33 with respect to p and setting ,6 = 0 
afterwards give 

(0)^ o = O - p + ^o) 2 (1 + 3 cos 2 e 0 ) , (42) 

P 0 in Equation 42 having the same meaning as that in Equation 41. 

(2) Let the distance in the x-direction from m 3 to the double point Lj be repre- 
sented by p . Then by a similar approximation, as before, Equation 35 reduces to 


(40) 

(41) 


v 1 “ LL U 


1 + ( 1 +3 cos 28 0 ) p 

By exactly the same procedure as before, we obtain 

(%) =E 1 (l + 3cos2 e 0 ), 

' H, 0=o 


-V = o. 


where 


El = 


(1 - p - p 0 ) p 0 


3 P 0 3 + 2m (1 + p 0 + p o 2 ) 

(^■) = (1 " Pn ~ P) 2 (1+3 cos 2 8 0 ) , 


(43) 


(44) 


(45) 

(46) 


and p 0 is the solution of Equation 43 with P = 0. 


I 
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(3) Let the distance in the x -direction from m 2 to the double point l 3 be repre- 
sented by l - p. Then Equation 35 reduces by approximation to 


[l + (1 + 3 cos 20 0 ) /5lJ (1 + n - p) - 1 1 _ - — _^ )2 ~ 0 • ( 47 ) 

From this is derived 


(If),.. = E 3<‘ * 3 '»* 


(48) 


where 


e 3 


(1 ■+ M - P 0 )(2 - P 0 ) 3 

3(2 - Po) 3 + 2 m [( 2 - Po> 2 + (2 - P 0 > + X J ’ 


(49) 


(3§) _ 0 = (! “ <°o + M) 2 (1 + 3 cos 20 o ) , (50) 

where p 0 is the solution of Equation 47 with /3 = 0 . 

Now if the positions and their corresponding values of C of the three double points 
on the x-axis are known for the case a = 0, their positions and the corresponding values 
of c for the case of small /? may be derived by 

• < 51 > 

and 

C = C " * ■ (“) 

The change in the ordinates of these points is given by Equation 38. 

For the earth-moon system, m = 0.01216. The values for the relevant quantities in 
this case are listed in Table 1. The derivative of the coordinates of these three points 

Table 1 

Values of p, c, e a , and Fj for the changes in 
position of L 1? L 2 , L 3 with 0 = 0 
[m = 0.01216] 



L x (i = 1) 

L 2 (i = 2) 

L 3 ( i = 3) 

p 

0.1510 

0.1679 

0.00709 

c 

3.18843 

3.17223 

3.01216 

E i 

0.0741 

-0.1566 

0.3326 

F i 

0.6053 

1.5831 

1 

-0.3820 



2Q 


3 


* 



(I) 


m 2 


55° 

0°^90°; 

L l 



m 


3 


( 2 ) 

N 

90 °— 0 ° 

55° m 2 


m 


3 M 


90° 

55° 


Figure 1 - Changes in position with /S at J3 - 0 of (1) the three double points, Lj , 

L 3 , and of (2) the inter secting points M and N with the x-axis of the zero-velocity surface 
passing through Lj. Notice the opposite directions of the changes in position of L 2 and M . 


with respect to p at their normal positions (s = 0) are furthermore illustrated in Fig- 
ure 1 for three positions of the sun (e Q - o, 77 / 4 , v/2 ). The directions and magnitudes 
of the arrows in the figure indicate the derivatives of the coordinates of these points with 
respect to p for three values of & c (0, v/4, v/2). 

In order to examine the change in position of the double points l 4 and l 5 , which 
make two equilateral triangles with m 2 and m 3 when p Q 0, we must resort to the 
original Equations 35 and 36. Differentiate them with respect to 5 , and set in the 
resulting equations: 


P - 0. r 2 = 

The required quantities (dx/d/3) 3=0 
simultaneously. 



and (dy/d/3) 3=0 are derived by solving the equations 


DEGENERATION OF THE CRITICAL 
ZERO-VELOCITY SURFACES 

When j3 = 0, which corresponds to the restricted three-body problem, the zero- 
velocity surface that passes through L 3 is frequently known as the inner contact surface 
and that which passes through L 2 , the outermost contact surface. The former intersects 
the x-axis at two more points besides Lj. Call the intersecting point on the positive x- 
axis M and that on the negative x-axis N. From the change in position of M and N with 
p, the general behavior of the system of zero-velocity surfaces can be inferred. 



10 


(1) Point M : Let its distance from m 3 be <j. Thus, from Equation 33, 

C = (l-/x + ^) 2 [l + (l + 3 cos 2^ 0 )/3] + ~ • 


(53) 


Differentiating Equation 53 with respect to 3 and utilizing the relation given by Equation 
46 give 


( 3?)». 0 = W‘ * 3 “■ “•> • 


(54) 


where 


E 


m 


2 


(p 0 + CT o)^ 2 - 2 m + °o _ Po) 


1 


- P + a o 


1 - M 

(1 + CT o ) 2 



(55) 


in which o-„ is the solution of Equation 53 with 3 = 0 and C = C 0 , corresponding to the 
inner contact surface of the restricted three-body problem. 


(2) Point n : Let its distance from m 2 be cr. Following the same procedure as 
before, we derive 



E n (l +3 cos 26 „) , 


where 


(1 - 2m - P 0 - + Po + ^o) 


M + a o 


1 - M 


(1 + °- 0 ) 2 J 


(56) 


(57) 


where ct 0 is cr of N when 3 = 0 . In both Equations 55 and 57, p 0 is the solution of Equa- 
tion 43 with 3 = 0. 

For M = 0.01216, 


E ra = 0.6217, and E n = 0.0589. (58) 

The changes in position of M and N are illustrated in Figure 1, from which it is seen 
that M moves out while L 2 moves in as 3(1 + 3 cos 2 6 0 ) increases. In other words, the 
inner contact surface will eventually meet the outermost contact surface at a certain 
value of 3(1 + 3 cos 26 0 ). When this happens, the two surfaces degenerate into one sur- 
face. It is evident from Equations 40, 54, 58 and Table 1 that the smallest value of 3 for 
which the critical surfaces become degenerated occurs at 9 0 = 0 . This threshold value 
of 3 (denoted by 3 C hereafter) can be determined in the following way: First calculate 
the two points L : and L 2 by Equation 35 with 6 0 = 0 . Denote the distances of I - 1 and L 2 
from m 3 by p x and p 2 , from which the corresponding values of C (denoted by Cj and 
C 2 , respectively) can be obtained from Equation 33. The degenerated case is given by 
the condition 
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Ci = c 2 , (59) 

which gives the required value /3 C . In Table 2 there is computed for the case p = 0.01216 
three sets of values from which we obtain 

/3 c = 0.0064 (60) 

by graphical interpolation. The corresponding values for p v p 2 = ^ M , and C x = C 2 are 
given in the last row of the table. Figure 2 illustrates the degenerated zero-velocity sur- 
face that passes through both L x and l 2 . 


Table 2 

Determination of 

[m = 0.01216] 


p 

Pi 

Ci 

?2 

C 2 

ct m 

0.0 

0.15097 

3.18843 

0.16788 

3.17223 

0.12580 

0.0025 

0.15171 

3.19542 

0.16672 

3.18557 

0.13306 

0.005 

0.15246 

3.20241 

0.16482 

3.19888 

0.14442 

0.0075 

0.15321 

3.20938 

0.16334 

3.21214 

(No solution) 

Pc = 

Pc = 

Pc = 

Pc - 

Pc = 

Pc = 

0.0064 

0.15288 

3.20632 

0.16398 

3.20632 

0.16398 


Figure 2 - Degenerate zero- 
velocity surface passing through 
both Lj and L 2 ; plotted for the 
case p =0.01216 corresponding 
to the earth-moon system 
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Further increase in /3 causes a fundamental change in shape of the critical zero- 
velocity surfaces. The surface that passes through L x opens up at M as can be seen 
from Equation 53, which does not yield any significant solution when > P c . The general 
behavior of the critical zero-velocity surfaces for the cases /3 > p c is illustrated in Fig- 
ure 3. In general the zero-velocity surfaces in the very restricted four-body problem 
are not symmetric even with respect to the x-axis. However, when o 0 = 0, they become 
symmetrical as is shown in both Figures 2 and 3. 



Figure 3 - Zero-velocity surfaces passing through 
hi and L 2 , respectively, when fi > fi c \ plotted for 
the case = 0.01216, P = 0.025, and d Q = 0 (fig- 
ure is symmetric with respect to the x-axis only 
because @ 0 = 0) 


DISCUSSION 

For the earth-moon-sun system, from Equation 34, 

P = 0.0028, (61) 

which is smaller than the threshold value for degeneracy as given by Equation 60. Thus, 
the critical surfaces will never become degenerated for any value of 6 0 . In other words, 
the inner and outermost contact surfaces for the earth-moon system can still be defined 
in spite of the presence of the sun, justifying the treatment in the previous paper of sat- 
ellites in the earth-moon-sun system as restricted three-body problems. 



For P < p c , the inner and outermost contact surfaces may be regarded as oscillating 
when 0 0 varies periodically. It follows from Equation 46 that a positive value of 
>6(1+3 cos 2& 0 ) makes a satellite escape easier than does a negative value. For example, 
a satellite with c > 3. 18843 will not escape from the neighborhood of the earth (or of the 
moon) in the framework of the restricted three-body problem (i.e., /3 = o). By the intro- 
duction of the fourth body (mj), a satellite will be retained inside the inner contact sur- 
face permanently only if c > 3. 19625. Similarly, the limiting value of c for retaining a 
satellite inside the outermost contact surface is now 3.18714, against 3.17223 in the 
restricted three-body problem. 

Although the present method of approach does not give the perturbation of orbital ele- 
ments of artificial satellites, it gives a general idea of where they could or could not go 
under given initial conditions when they are no longer very near to the earth. 


NASA - Langley Field, Va. 



